Hybrid Finite Volume Scheme for a Two-phase Flow in Heterogeneous Porous Media

We propose a finite volume method on general meshes for the numerical simulation of an incompressible and immiscible two-phase flow in porous media. We consider the case that can be written as a coupled system involving a degenerate parabolic convection-diffusion equation for the saturation together with a uniformly elliptic equation for the global pressure. The numerical scheme, which is implicit in time, allows computations in the case of a heterogeneous and anisotropic permeability tensor. The convective fluxes, which are non monotone with respect to the unknown saturation and discontinuous with respect to the space variables, are discretized by means of a special Godunov scheme. We prove the existence of a discrete solution which converges, along a subsequence, to a solution of the continuous problem. We present a number of numerical results in space dimension two, which confirm the efficiency of the numerical method. Résumé. Nous proposons un schéma de volumes finis hybrides pour la discrétisation d’un problème d’écoulement diphasique incompressible et immiscible en milieu poreux. On suppose que ce problème a la forme d’une équation parabolique dégénérée de convection-diffusion en saturation couplée à une équation uniformément elliptique en pression. On considère un schéma implicite en temps, où les flux diffusifs sont discrétisés par la méthode des volumes finis hybride, ce qui permet de pouvoir traiter le cas d’un tenseur de perméabilité anisotrope et hétérogène sur un maillage très général, et l’on s’appuie sur un schéma de Godunov pour la discrétisation des flux convectifs, qui peuvent être non monotones et discontinus par rapport aux variables spatiales. On démontre l’existence d’une solution discrète, dont une sous-suite converge vers une solution faible du problème continu. On présente finalement des cas test bidimensionnels. Introduction We study the simplified dead-oil model in so-called global pressure or fractional flow formulation [5]. Let Ω be a polyhedral open bounded connected subset of R and T > 0; we consider the following system u = −K (λ(s)∇p− ξ(s)g) , −∇ · u = qw + qo in Ω× (0, T ), (1) ω ∂s ∂t +∇ · (uf(s) + γ(s)Kg)−∇ · (K∇φ(s)) = qo in Ω× (0, T ), (2) together with some initial condition and homogeneous Dirichlet boundary conditions s(·, 0) = s0 in Ω, p = 0 and s = 0 on ∂Ω× (0, T ). (3) ∗ This work was supported by the GdR MoMaS (PACEN/CNRS, ANDRA, BRGM, CEA, EdF, IRSN), France. 1 Laboratoire de Mathématiques, Université de Paris-Sud 11, F-91405 Orsay Cedex, France c © EDP Sciences, SMAI 2012 Article published online by EDP Sciences and available at http://www.esaim-proc.org or http://dx.doi.org/10.1051/proc/201235016 ESAIM: PROCEEDINGS 211 The unknown functions are the global pressure p and the saturation of the oil-phase s. The usual assumption λo + λw ≥ λ > 0 implies that the first equation is uniformly elliptic in p, whereas the second one is parabolic degenerate with respect to s. The heterogeneity and anisotropy of porous media is a numerical challenge even when studying the elliptic problem derived from Darcy’s law for one-phase problem. Many schemes were proposed and analyzed in the last decades for its discretization. See [9] and [7] for more references and for the detailed description and comparison of those numerical methods. In this paper we propose an implicit fully coupled hybrid finite volume scheme for the problem (1)-(3). The discretization of the diffusion terms is based upon a hybrid finite volume method [6], which allows the tensor K to be anisotropic and highly variable in space. Remark that heterogeneity also affects the first order terms K(x)gξ(·) and uf(·) +K(x)gγ(·), which may be discontinuous with respect to the space variable x; in that case they require a suitable treatment see e.g. [8, 12, 13]. We apply the Godunov scheme proposed in [10], which seems to be a natural choice, since the hybrid (interface) unknowns are used. Assumptions on the data: (H1) φ ∈ C(R), φ(0) = 0, is a strictly increasing piecewise continuously differentiable Lipschitz continuous function. We assume that the function φ is Hölder continuous, namely that there exists Hφ > 0 and α ∈ (0, 1] such that |s1− s2| ≤ Hφ|φ(s1)−φ(s2)|. It is also such that φ is Lipschitz continuous on R \ (0, 1); (H2) The functions λ, ξ, γ, f ∈ C(R) are Lipschitz continuous; (H2a) ξ and γ are convex functions, moreover γ is such that γ(0) = γ(1) = 0; g is a constant vector from R; (H2b) f is a nondecreasing function and it satisfies f(0) = 0, f(1) = 1; (H3) s0 ∈ L(Ω); ω ∈ L(Ω); (H4) qo, qw ∈ L(0, T ;L(Ω)) and such that qo + qw ≥ 0 a.e. in Ω× (0, T ). The physical range of values for the saturation is [0, 1], however we extend the definition of all nonlinear functions in (1)-(2) outside of [0, 1], since the numerical scheme which we study does not preserve neither the maximum principle, nor the positivity of s (nor the bound s ≤ 1). Assumptions on the geometry: (H5a) Ω is a polyhedral open bounded connected subset of R, with d ∈ N, and ∂Ω = Ω\Ω its boundary. (H5b) K is a piecewise constant function from Ω to Md(R), where Md(R) denotes the set of real d×d matrices. More precisely we assume that there exist a finite family (Ωi)i∈{1,...,I} of open connected polyhedral in R , such that Ω = ⋃ i∈{1,...,I}Ωi, Ωi ⋂ Ωj = ∅ if i 6= j and such that K(x)|Ωi = Ki ∈ Md(R). By Γi,j we denote the interface between the sub-domains i, j, Γi,j = Ωi ⋂ Ωj . We assume that there exist two positive constants K and K such that the eigenvalues of the symmetric positive definite Ki are included in [K,K] for all i ∈ {1, . . . , I}. For the sake of simplicity we have assumed that the heterogeneity of the medium is only expressed through the x-dependence of the absolute permeability tensor K = K(x). However it is very simple to extend the analysis to the case where λ, ξ, γ and f depend on the rock type. On the other hand if we suppose that φ is discontinuous in space, this may lead to significant difficulties. The analysis of the case where the capillary pressure field (and also φ) is discontinuous was carried out in [4] and [3]. It is also worth noting that the partitioning of Ω introduced in H5b is only used in order to provide a control on the gravity terms. In the case that the gravity effects are neglected, one can consider a fully heterogenous permeability field K. 1. The finite volume scheme 1.1. The main definitions Definition 1.1 (Discretization of Ω). Let Ω be a polyhedral open bounded connected subset of R, with d ∈ N, ∂Ω = Ω\Ω its boundary, and (Ωi)i∈{1,...,I} is the partition of Ω in the sense of (H5b). A discretization of Ω, denoted by D, is defined as the triplet D = (M, E ,P), where: 1. M is a finite family of non empty connex open disjoint subsets of Ω (the ”control volumes”) such that Ω = ⋃ K∈M K. For any K ∈ M, let ∂K = K\K be the boundary of K; we define m(K) > 0 as the measure of K and hK as the diameter of K. We also assume that the mesh resolve the structure of the medium, i.e. 212 ESAIM: PROCEEDINGS for all K ∈ M there exist i ∈ {1, . . . , I} such that K ⊂ Ωi. The size of the discretization D is defined by hD = supK∈M diam(K). 2. E is a finite family of disjoint subsets of Ω (the ”edges” of the mesh), such that, for all σ ∈ E, σ is a non empty open subset of a hyperplane of R, whose (d− 1)-dimensional measure m(σ) is strictly positive. We also assume that, for all K ∈ M, there exists a subset EK of E such that ∂K = ⋃ σ∈EK σ. For each σ ∈ E, we set Mσ = {K ∈ M|σ ∈ EK}. We then assume that, for all σ ∈ E, either Mσ has exactly one element and then σ ∈ ∂Ω (the set of these interfaces called boundary interfaces, is denoted by Eext) or Mσ has exactly two elements (the set of these interfaces called interior interfaces, is denoted by Eint). For all σ ∈ E, we denote by xσ the barycenter of σ. For all K ∈ M and σ ∈ EK, we denote by nK,σ the outward normal unit vector. 3. P is a family of points of Ω indexed by M, denoted by P = (xK)K∈M, such that for all K ∈ M, xK ∈ K; moreover K is assumed to be xK-star-shaped, which means that for all x ∈ K, there holds [xK ,x] ∈ K. Denoting by dK,σ the Euclidean distance between xK and the hyperplane containing σ, one assumes that dK,σ > 0. We denote by DK,σ the cone of vertex xK and basis σ. Definition 1.2 (Time discretization). We divide the time interval (0, T ) into N equal time steps of length δt = T/N , where δt is the uniform time step defined by δt = tn − tn−1. Definition 1.3 (The hybrid space XD(Ω)). Let D = (M, E ,P) be a discretization of Ω. We define XD = {((vK)K∈M, (vσ)σ∈E), vK ∈ R, vσ ∈ R}, XD,0 = {v ∈ XD such that (vσ)σ∈Eext = 0}. (4) Taking into account the time discretization we define XD,δt = X N D and XD,δt,0 = X N D,0. 1.2. The numerical scheme Let us introduce the discrete saturation ( (sK)K∈M , (s n σ)σ∈E ) n∈{1,...,N} ∈ XD,δt and the discrete global pressure ( (pK)K∈M , (p n σ)σ∈E ) n∈{1,...,N} ∈ XD,δt, which are the main discrete unknowns. Let f denote λ, ξ, γ or f , we introduce the following notation fK = f(s n K) and f n σ = f(s n σ). Next, let q n i,K denote the mean value of the source terms, i.e. q i,K = 1 m(K)δt ∫ tn tn−1 ∫ K qi(x, t) dxdt with i ∈ {w, n} and ω(K) denotes the porous volume of the element K, ω(K) = ∫ K ω(x) dx. Let U K,σ be an approximation of the total flux through the interface σ U K,σ ≈ 1 δt ∫ tn tn−1 ∫ σ K(λ(s)∇p − ξ(s)g) · nK,σdνdt (5) and let F K,σ be an approximation of the non-wetting phase flux F K,σ ≈ 1 δt ∫ tn tn−1 ∫ σ (uf(s) + γ(s)Kg −K∇φ(s)) · nK,σdνdt. (6) The numerical fluxes U K,σ and F n K,σ have to be constructed as functions of the discrete unknowns. Following the ideas of [6] we write the scheme in the variational form. For each n ∈ {1, . . . , N} find s ∈ XD,0 and p ∈ XD,0 such that for all v, w ∈ XD,0: